home *** CD-ROM | disk | FTP | other *** search
/ HAKERIS 11 / HAKERIS 11.ISO / linux / system / LinuxConsole 0.4 / linuxconsole0.4install-en.iso / guile0.4.lcm / share / guile / slib / root.scm < prev    next >
Encoding:
Text File  |  2004-01-06  |  7.1 KB  |  218 lines
  1. ;;;"root.scm" Newton's and Laguerre's methods for finding roots.
  2. ;Copyright (C) 1996, 1997 Aubrey Jaffer
  3. ;
  4. ;Permission to copy this software, to redistribute it, and to use it
  5. ;for any purpose is granted, subject to the following restrictions and
  6. ;understandings.
  7. ;
  8. ;1.  Any copy made of this software must include this copyright notice
  9. ;in full.
  10. ;
  11. ;2.  I have made no warrantee or representation that the operation of
  12. ;this software will be error-free, and I am under no obligation to
  13. ;provide any services, by way of maintenance, update, or otherwise.
  14. ;
  15. ;3.  In conjunction with products arising from the use of this
  16. ;material, there shall be no use of my name in any advertising,
  17. ;promotional, or sales literature without prior written consent in
  18. ;each case.
  19.  
  20. (require 'logical)
  21.  
  22. ;;;; Newton's Method explained in:
  23. ;;; D. E. Knuth, "The Art of Computer Programming", Vol 2 /
  24. ;;; Seminumerical Algorithms, Reading Massachusetts, Addison-Wesley
  25. ;;; Publishing Company, 2nd Edition, p. 510
  26.  
  27. (define (newton:find-integer-root f df/dx x_0)
  28.   (let loop ((x x_0) (fx (f x_0)))
  29.     (cond
  30.      ((zero? fx) x)
  31.      (else
  32.       (let ((df (df/dx x)))
  33.     (cond
  34.      ((zero? df) #f)        ; stuck at local min/max
  35.      (else
  36.       (let* ((delta (quotient (+ fx (quotient df 2)) df))
  37.          (next-x (cond ((not (zero? delta)) (- x delta))
  38.                    ((positive? fx) (- x 1))
  39.                    (else (- x -1))))
  40.          (next-fx (f next-x)))
  41.         (cond ((>= (abs next-fx) (abs fx)) x)
  42.           (else (loop next-x next-fx)))))))))))
  43.  
  44. (define (integer-sqrt y)
  45.   (newton:find-integer-root (lambda (x) (- (* x x) y))
  46.                 (lambda (x) (* 2 x))
  47.                 (ash 1 (quotient (integer-length y) 2))))
  48.  
  49. (define (newton:find-root f df/dx x_0 prec)
  50.   (if (and (negative? prec) (integer? prec))
  51.       (let loop ((x x_0) (fx (f x_0)) (count prec))
  52.     (cond ((zero? count) x)
  53.           (else (let ((df (df/dx x)))
  54.               (cond ((zero? df) #f) ; stuck at local min/max
  55.                 (else (let* ((next-x (- x (/ fx df)))
  56.                      (next-fx (f next-x)))
  57.                     (cond ((= next-x x) x)
  58.                       ((> (abs next-fx) (abs fx)) #f)
  59.                       (else (loop next-x next-fx
  60.                               (+ 1 count)))))))))))
  61.       (let loop ((x x_0) (fx (f x_0)))
  62.     (cond ((< (abs fx) prec) x)
  63.           (else (let ((df (df/dx x)))
  64.               (cond ((zero? df) #f) ; stuck at local min/max
  65.                 (else (let* ((next-x (- x (/ fx df)))
  66.                      (next-fx (f next-x)))
  67.                     (cond ((= next-x x) x)
  68.                       ((> (abs next-fx) (abs fx)) #f)
  69.                       (else (loop next-x next-fx))))))))))))
  70.  
  71. ;;; H. J. Orchard, "The Laguerre Method for Finding the Zeros of
  72. ;;; Polynomials", IEEE Transactions on Circuits and Systems, Vol. 36,
  73. ;;; No. 11, November 1989, pp 1377-1381.
  74.  
  75. (define (laguerre:find-root f df/dz ddf/dz^2 z_0 prec)
  76.   (if (and (negative? prec) (integer? prec))
  77.       (let loop ((z z_0) (fz (f z_0)) (count prec))
  78.     (cond ((zero? count) z)
  79.           (else
  80.            (let* ((df (df/dz z))
  81.               (ddf (ddf/dz^2 z))
  82.               (disc (sqrt (- (* df df) (* fz ddf)))))
  83.          (if (zero? disc)
  84.              #f
  85.              (let* ((next-z
  86.                  (- z (/ fz (if (negative? (+ (* (real-part df)
  87.                                  (real-part disc))
  88.                               (* (imag-part df)
  89.                                  (imag-part disc))))
  90.                         (- disc) disc))))
  91.                 (next-fz (f next-z)))
  92.                (cond ((>= (magnitude next-fz) (magnitude fz)) z)
  93.                  (else (loop next-z next-fz (+ 1 count))))))))))
  94.       (let loop ((z z_0) (fz (f z_0)) (delta-z #f))
  95.     (cond ((< (magnitude fz) prec) z)
  96.           (else
  97.            (let* ((df (df/dz z))
  98.               (ddf (ddf/dz^2 z))
  99.               (disc (sqrt (- (* df df) (* fz ddf)))))
  100.          ;;(print 'disc disc)
  101.          (if (zero? disc)
  102.              #f
  103.              (let* ((next-z
  104.                  (- z (/ fz (if (negative? (+ (* (real-part df)
  105.                                  (real-part disc))
  106.                               (* (imag-part df)
  107.                                  (imag-part disc))))
  108.                         (- disc) disc))))
  109.                 (next-delta-z (magnitude (- next-z z))))
  110.                ;;(print 'next-z next-z )
  111.                ;;(print '(f next-z) (f next-z))
  112.                ;;(print 'delta-z delta-z 'next-delta-z next-delta-z)
  113.                (cond ((zero? next-delta-z) z)
  114.                  ((and delta-z (>= next-delta-z delta-z)) z)
  115.                  (else
  116.                   (loop next-z (f next-z) next-delta-z)))))))))))
  117.  
  118. (define (laguerre:find-polynomial-root deg f df/dz ddf/dz^2 z_0 prec)
  119.   (if (and (negative? prec) (integer? prec))
  120.       (let loop ((z z_0) (fz (f z_0)) (count prec))
  121.     (cond ((zero? count) z)
  122.           (else
  123.            (let* ((df (df/dz z))
  124.               (ddf (ddf/dz^2 z))
  125.               (tmp (* (+ deg -1) df))
  126.               (sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf))))
  127.               (df+sqrt-H (+ df sqrt-H))
  128.               (df-sqrt-H (- df sqrt-H))
  129.               (next-z
  130.                (- z (/ (* deg fz)
  131.                    (if (>= (magnitude df+sqrt-H)
  132.                        (magnitude df-sqrt-H))
  133.                    df+sqrt-H
  134.                    df-sqrt-H)))))
  135.          (loop next-z (f next-z) (+ 1 count))))))
  136.       (let loop ((z z_0) (fz (f z_0)))
  137.     (cond ((< (magnitude fz) prec) z)
  138.           (else
  139.            (let* ((df (df/dz z))
  140.               (ddf (ddf/dz^2 z))
  141.               (tmp (* (+ deg -1) df))
  142.               (sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf))))
  143.               (df+sqrt-H (+ df sqrt-H))
  144.               (df-sqrt-H (- df sqrt-H))
  145.               (next-z
  146.                (- z (/ (* deg fz)
  147.                    (if (>= (magnitude df+sqrt-H)
  148.                        (magnitude df-sqrt-H))
  149.                    df+sqrt-H
  150.                    df-sqrt-H)))))
  151.          (loop next-z (f next-z))))))))
  152.  
  153. (define (secant:find-root-1 f x0 x1 prec must-bracket?)
  154.   (letrec ((stop?
  155.         (cond ((procedure? prec) prec)
  156.           ((and (integer? prec) (negative? prec))
  157.            (lambda (x0 x1 fmax count)
  158.              (>= count (- prec))))
  159.           (else
  160.            (lambda (x0 f0 x1 f1 count)
  161.              (and (< (abs f0) prec)
  162.               (< (abs f1) prec))))))
  163.        (bracket-iter
  164.         (lambda (xlo flo glo xhi fhi ghi count)
  165.           (define (step xnew fnew)
  166.         (cond ((or (= xnew xlo)
  167.                (= xnew xhi))
  168.                (let ((xmid (+ xlo (* 1/2 (- xhi xlo)))))
  169.              (if (= xnew xmid)
  170.                  xmid
  171.                  (step xmid (f xmid)))))
  172.               ((positive? fnew)
  173.                (bracket-iter xlo flo (if glo (* 0.5 glo) 1)
  174.                      xnew fnew #f
  175.                      (+ count 1)))
  176.               (else
  177.                (bracket-iter xnew fnew #f
  178.                      xhi fhi (if ghi (* 0.5 ghi) 1)
  179.                      (+ count 1)))))
  180.           (if (stop? xlo flo xhi fhi count)
  181.           (if (> (abs flo) (abs fhi)) xhi xlo)
  182.           (let* ((fflo (if glo (* glo flo) flo))
  183.              (ffhi (if ghi (* ghi fhi) fhi))
  184.              (del (- (/ fflo (- ffhi fflo))))
  185.              (xnew (+ xlo (* del (- xhi xlo))))
  186.              (fnew (f xnew)))
  187.             (step xnew fnew))))))
  188.     (let ((f0 (f x0))
  189.       (f1 (f x1)))
  190.       (cond ((<= f0 0 f1)
  191.          (bracket-iter x0 f0 #f x1 f1 #f 0))
  192.         ((<= f1 0 f0)
  193.          (bracket-iter x1 f1 #f x0 f0 #f 0))
  194.         (must-bracket? #f)
  195.         (else
  196.          (let secant-iter ((x0 x0)
  197.                    (f0 f0)
  198.                    (x1 x1)
  199.                    (f1 f1)
  200.                    (count 0))
  201.            (cond ((stop? x0 f0 x1 f1 count)
  202.               (if (> (abs f0) (abs f1)) x1 x0))
  203.              ((<= f0 0 f1)
  204.               (bracket-iter x0 f0 #f x1 f1 #f count))
  205.              ((>= f0 0 f1)
  206.               (bracket-iter x1 f1 #f x0 f0 #f count))
  207.              ((= f0 f1) #f)
  208.              (else
  209.               (let* ((xnew (+ x0 (* (- (/ f0 (- f1 f0))) (- x1 x0))))
  210.                  (fnew (f xnew))
  211.                  (fmax (max (abs f1) (abs fnew))))
  212.             (secant-iter x1 f1 xnew fnew (+ count 1)))))))))))
  213.  
  214. (define (secant:find-root f x0 x1 prec)
  215.   (secant:find-root-1 f x0 x1 prec #f))
  216. (define (secant:find-bracketed-root f x0 x1 prec)
  217.   (secant:find-root-1 f x0 x1 prec #t))
  218.